###########
# STEP 2: EXPERIMENT ON SOCIAL VALUE OF PARTICIPATION
###########
# ---- Produces 
	# Figure 1
	# Figure 2
	# Appendix Table A1
	# Appendix Figure A1

###########
# STEP 3: SURVEY PLUS NEIGHBORHOOD DATA
###########
# ---- Produces 
	# Table 3
	# Figure 3


#-------Load Packages------#
#install.packages()
require(foreign)
require(xtable)
require(car)
require(lattice)
require(MASS)
require(plotrix)
require(survey)
require(gmodels)
require(prettyR)
require(RColorBrewer)
require(gplots) #For creating histograms with confidence intervals - barplot2
require(ggplot2)
require(plyr) #for ddply command
require("boot") #for creating weighted ses
require(stringr) #mangaing strings
require(weights)  #for wtd.t.test



###########
# ------- LOAD AND MERGE SURVEY DATA
###########


#----- Set working directory 
setwd("~/Box Sync/S_CollectiveActionPaper/Chapter2/PaperDrafts/Submissions/APSR_3_FinalSubmit/Replication/Data/SlimData/")


# -------- Load in three data sets: National Survey, Black Sample, Latino Sample
ns <- read.csv(file="2National_Slim.csv")  # National Survey Sample
bs <- read.csv(file="2Black_Slim.csv")  # Black Only Survey Sample
ls <- read.csv(file="2Latino_Slim.csv") # Latino Only Survey Sample


#------ Merge Data Sets into single data set

# Latino + National Sample: Merge and keep non-matching columns
data1 <- rbind.fill(ns, ls)
names(data1)
length(data1$caseid)
length(ns$caseid)
length(ls$caseid)

#Merge Latino/Nation + Black.  Keep non-matching columns. 
data2 <- rbind.fill(data1, bs)
length(data1$caseid)
length(data2$caseid)



###########
# ------- RECODING FUN
###########
data <- data2

#Recode Treatment Assignments: Community Actions
table(data$apa_comact1)
data$apa_comact1<-str_trim(data$apa_comact1) #strips white space from end of variable name to match categories
table(data$apa_comact1)
data$p1act <- recode(data$apa_comact1, "
					'Annually attends charitable events' = 1;
					'Annually attends political rallies' = 2;
					'Votes in every presidential election' = 3
					")
table(data$p1act)

table(data$apa_comact2)
data$apa_comact2<-str_trim(data$apa_comact2)
table(data$apa_comact2)
data$p2act <- recode(data$apa_comact2, "
					'Annually attends charitable events' = 1;
					'Annually attends political rallies' = 2;
					'Votes in every presidential election' = 3
					")
table(data$p2act)

data$apa_comact3<-str_trim(data$apa_comact3)
table(data$apa_comact3)
data$p3act <- recode(data$apa_comact3, "
					'Annually attends charitable events' = 1;
					'Annually attends political rallies' = 2;
					'Votes in every presidential election' = 3
					")
table(data$p3act)


# Recode Dependent Variables

#----------Likability
#----------Vote
data$vlike<- rep(NA, length(data$caseid))
data$vlike[data$p2act==3]<- data$APA002[data$p2act==3]
data$vlike[data$p1act==3]<- data$APA004[data$p1act==3]
data$vlike[data$p3act==3]<- data$APA007[data$p3act==3]
table(data$vlike)

#Rally
data$rlike <- rep(NA, length(data$caseid))
data$rlike[data$p2act==2]  <- data$APA002[data$p2act==2] 
data$rlike[data$p1act==2] <- data$APA004[data$p1act==2] 
data$rlike[data$p3act==2] <- data$APA007[data$p3act==2] 
table(data$rlike)

#Charity
data$clike <- rep(NA,  length(data$caseid))
data$clike[data$p2act==1]  <- data$APA002[data$p2act==1] 
data$clike[data$p1act==1] <- data$APA004[data$p1act==1] 
data$clike[data$p3act==1] <- data$APA007[data$p3act==1] 
table(data$clike)

#---Baseline/Potluck
data$blike <- data$APA010



#------- RESPECTABILITY SCORE
#Vote
data$v.resp <- rep(NA, length(data$caseid))
data$v.resp[data$p2act==3]  <- data$APA003[data$p2act==3] 
data$v.resp[data$p1act==3]  <- data$APA005[data$p1act==3] 
data$v.resp[data$p3act==3]  <- data$APA008[data$p3act==3] 
table(data$v.resp)

#Rally
data$r.resp <- rep(NA, length(data$caseid))
data$r.resp[data$p2act==2]  <- data$APA003[data$p2act==2] 
data$r.resp[data$p1act==2]  <- data$APA005[data$p1act==2] 
data$r.resp[data$p3act==2]  <- data$APA008[data$p3act==2] 
table(data$r.resp)

#Charity
data$c.resp <- rep(NA, length(data$caseid))
data$c.resp[data$p2act==1] <- data$APA003[data$p2act==1] 
data$c.resp[data$p1act==1] <- data$APA005[data$p1act==1] 
data$c.resp[data$p3act==1] <- data$APA008[data$p3act==1] 
table(data$c.resp)

#Baseline/Potuck
data$b.resp <- data$APA011
table(data$b.resp)


#------Recode Race
names(data)
table(data$race)
data$race_sm<- rep(NA, length(data$caseid))
data$race_sm[data$race=="White"] <- "White"
data$race_sm[data$race=="Black"] <- "Black"
data$race_sm[data$race=="Hispanic"] <- "Hispanic"
data$race_sm[data$race=="Asian"] <- "Asian"

#Race without Asians
table(data$race_sm)
data$race3 <- rep(NA, length(data$caseid))
data$race3[data$race_sm=="White"] <-"White"
data$race3[data$race_sm=="Black"] <-"Black"
data$race3[data$race_sm=="Hispanic"] <-"Hispanic"
table(data$race3)

#Race to match census categories
table(data$race)
data$demRace <- recode(data$race, " 'White' = 1; 
						'Black' = 2; 
						'Native American' = 3; 
						'Asian' = 4; 
						'Middle Eastern' = 5;
						'Mixed' = 5;
						'Other' = 5;
						else = NA   ")
table(data$demRace)

#Hispanic
table(data$hisp)
data$demHisp<- rep(2, length(data$caseid)) #2 indicates not hispanic
data$demHisp[data$race=="Hispanic"] <- 1 #combine hispanic from three sources
data$demHisp[data$hisp == "Yes"] <- 1
table(data$demHisp)

#Race and Hispanic Combined
table(data$demRace, data$demHisp)
data$demRaceFull <-  rep(NA, length(data$caseid))
data$demRaceFull[data$demRace == 1] <- 1 #White, not hispanic
data$demRaceFull[data$demRace == 2 & data$demHisp == 2 ] <- 2 #Black, not hispanic
data$demRaceFull[data$demRace == 3] <- 3 #Native American, not hispanic
data$demRaceFull[data$demRace == 4] <- 4 #Asian, not hispanic
data$demRaceFull[data$demRace == 5] <- 5 #Other, not hispanic
data$demRaceFull[data$demHisp == 1] <- 6 #Hispanic, any race
table(data$demRaceFull)
table(data$demHisp)

#Gender
table(data$gender)
data$demGender <-  recode(data$gender, " 'Male' = 1; 'Female' = 2" )
table(data$demGender)

data$edu_sm <- rep(NA, length(data$caseid))
data$edu_sm[as.numeric(data$edu1)==1] <-1
data$edu_sm[as.numeric(data$edu1)==2] <-1
data$edu_sm[as.numeric(data$edu1)==3] <-2
data$edu_sm[as.numeric(data$edu1)==4] <-2
data$edu_sm[as.numeric(data$edu1)==5] <-3
data$edu_sm[as.numeric(data$edu1)==6] <-3
table(data$edu_sm)

#Into Four Components
data$edu2 <- rep(NA, length(data$caseid))
data$edu2[as.numeric(data$edu1)==1] <-1
data$edu2[as.numeric(data$edu1)==2] <-2
data$edu2[as.numeric(data$edu1)==3] <-2
data$edu2[as.numeric(data$edu1)==4] <-2
data$edu2[as.numeric(data$edu1)==5] <-3
data$edu2[as.numeric(data$edu1)==6] <-4
table(data$edu2)

#Matching CPS weight categories; #recode to 5 categories
table(data$edu1)
data$edu5 <- recode(data$edu1, " 1 = 1; 2 = 2;  3= 3; 4 = 4; 5 = 4; 6 = 5   ")
table(data$edu5)



#Qunatile test by thirds
quantile(as.numeric(data$inc), c(0.3333333, 0.6666666), na.rm=TRUE)
data$inc_sm <- rep(NA, length(data$caseid))
data$inc_sm[as.numeric(data$inc)==1]<-1
data$inc_sm[as.numeric(data$inc)==2]<-1
data$inc_sm[as.numeric(data$inc)==3]<-1
data$inc_sm[as.numeric(data$inc)==4]<-2
data$inc_sm[as.numeric(data$inc)==5]<-2
data$inc_sm[as.numeric(data$inc)==6]<-2
data$inc_sm[as.numeric(data$inc)==7]<-3
data$inc_sm[as.numeric(data$inc)==8]<-3
data$inc_sm[as.numeric(data$inc)==9]<-3
data$inc_sm[as.numeric(data$inc)==10]<-3
data$inc_sm[as.numeric(data$inc)==11]<-3
data$inc_sm[as.numeric(data$inc)==12]<-3
#Low == 1 (0-29,999)
#Medium == 2 (30,000 - 59,99)
#High == 3 (60,000 +)
table(data$inc_sm)
sum(table(data$inc_sm))

#Matching income categories to overperformance paper
data$inc2 <- rep(NA, length(data$caseid))
data$inc2[as.numeric(data$inc)==1]<-1
data$inc2[as.numeric(data$inc)==2]<-1
data$inc2[as.numeric(data$inc)==3]<-2
data$inc2[as.numeric(data$inc)==4]<-2
data$inc2[as.numeric(data$inc)==5]<-3
data$inc2[as.numeric(data$inc)==6]<-3
data$inc2[as.numeric(data$inc)==7]<-3
data$inc2[as.numeric(data$inc)==8]<-3
data$inc2[as.numeric(data$inc)==9]<-4
data$inc2[as.numeric(data$inc)==10]<-4
data$inc2[as.numeric(data$inc)==11]<-4
data$inc2[as.numeric(data$inc)==12]<-4
table(data$inc2)

#Recode Income to Mach CPS file for weights
table(data$inc)
data$demInc <- rep(NA, length(data$caseid))
data$demInc[data$inc == 1] <-1
data$demInc[data$inc == 2] <-  1
data$demInc[data$inc == 3] <- 2
data$demInc[data$inc == 4] <-2
data$demInc[data$inc == 5] <-  2
data$demInc[data$inc == 6] <- 3
data$demInc[data$inc == 7] <- 3
  data$demInc[data$inc == 8] <- 3
data$demInc[data$inc == 9] <- 4
data$demInc[data$inc == 10] <- 5
data$demInc[data$inc == 11] <- 5
data$demInc[data$inc == 12] <- 5
table(data$demInc)

#Make combined ses measure
data$ses <- rep(NA, length(data$caseid))
data$ses[data$inc_sm ==1 & data$edu_sm ==1]<-1
data$ses[data$inc_sm ==2 & data$edu_sm ==2]<-2
data$ses[data$inc_sm ==3 & data$edu_sm ==3]<-3
table(data$ses)
sum(table(data$ses)) #looses half of cases because there are some folks who are low ses high edu etc

#Recode party id
table(data$pid3)
data$pid_clean <- rep(NA, length(data$caseid))
data$pid_clean[data$pid3 == "Democrat" ]<- "dem"
data$pid_clean[data$pid3 == "Republican" ]<- "rep"
data$pid_clean[data$pid3 == "Independent" ]<- "ind"
table(data$pid_clean)

#Age
#Population is over 18
table(data$birthyr)
data$age <- 2015 - data$birthyr
table(data$age)

#Make Age categories Match CPS
data$ageDem <- rep(NA, length(data$caseid))
data$ageDem[data$age  >18 & data$age < 25  ] <- 	1	#  Between 18 and 24 - {Between 18 and 24}
data$ageDem[data$age  > 24 & data$age < 35 ] <- 	2	#  Between 25 and 34 - {Between 25 and 34}
data$ageDem[data$age  > 34 & data$age < 45 ] <- 	3#  Between 35 and 44 - {Between 35 and 44}
data$ageDem[data$age  > 44 & data$age < 55 ] <- 	4#  Between 45 and 54 - {Between 45 and 54}
data$ageDem[data$age  > 54 & data$age < 65 ] <- 	5#  Between 55 and 64 - {Between 55 and 64}
data$ageDem[data$age  > 64 & data$age < 75 ] <- 	6#  Between 65 and 74 - {Between 65 and 74}
data$ageDem[data$age  > 74  ] <- 	7#  75 +
table(data$ageDem)


#Region 
#Merge in FIPS State code and regiona division information
statecode<- read.csv(file="2FIPS_statecode.csv")
names(statecode)
statecode$demState <- statecode$fips_state
statecode$state2 <- statecode$name
statecode$div9 <- statecode$census_division
table(statecode$div9, statecode$census_division_name)


#merge in state codes and regions
data3 <- merge(data, statecode, by = "state2") 
names(data3)






###########
# ------- CONSTRUCT WEIGHT
###########

#Remove observations missing income in preparation for weight construction 
table(is.na(data3$demInc))
data4 <- data3[complete.cases(data3[, c("demInc")]), ]
length(data4$caseid)


## Read in CPS Data ----------
cps <- suppressWarnings(read.spss("2CPS_Weight.sav", to.data.frame = TRUE, use.value.labels = FALSE))


## Format Data variables --------------

#Population Weights
cps$wts_cps <- cps$PWCMPWGT  # Final composite weight, for manufacturing 
#Race
cps$demRace <- cps$RECODE2 #(1 = White, 2 = Black, 3 = American Indian, 4= Asian, 5 = Other)
table(cps$demRace)

#Hispanic
cps$demHisp <- cps$PEHSPNON

#Create combined race and hispanic measure, like in origonal data set
cps$demRaceFull <- rep(NA, length(cps$GESTFIPS))
cps$demRaceFull[cps$demRace == 1 & cps$demHisp == 2 ] <- 1 #White, not hispanic
cps$demRaceFull[cps$demRace == 2 & cps$demHisp == 2 ] <- 2 #Black, not hispanic
cps$demRaceFull[cps$demRace == 3 & cps$demHisp == 2 ] <- 3 #Native American, not hispanic
cps$demRaceFull[cps$demRace == 4 & cps$demHisp == 2 ] <- 4 #Asian, not hispanic
cps$demRaceFull[cps$demRace == 5 & cps$demHisp == 2 ] <- 5 #Other, not hispanic
cps$demRaceFull[cps$demHisp == 1 ] <- 6 #Hispanic, any race
table(cps$demRaceFull)

#Gender
cps$demGender <- cps$PESEX
table(cps$demGender)

#Age
cps<-  subset(cps, cps$RECODE1 != 1)
table(cps$RECODE1)
cps$ageDem <- recode(cps$RECODE1, "2=1;3=2;4=3;5=4;6=5;7=6;8=7;else=NA")
table(cps$ageDem)
levels(cps$ageDem) <- levels(cps$RECODE1)[2:8]

#Education
#what do I do with the 6th category here!??
cps$edu5 <- recode(cps$RECODE3, "6 = NA")
table(cps$edu5)
levels(cps$edu5) <- levels(cps$RECODE3)[1:5]

#Income
cps$demInc <- cps$RECODE4
table(cps$demInc)

#Region
cps$demState <- as.numeric(cps$GESTFIPS)
head(cps$demState)
cps$div9 <- rep(NA, dim(cps)[1]) #four regions - northeast, midwest, south, west
cps$div9[cps$demState %in% c(20, 30, 46, 22, 40, 7)] <- 1
cps$div9[cps$demState %in% c(33, 39, 31)] <- 2
cps$div9[cps$demState %in% c(50, 23, 14, 15, 36)] <- 3
cps$div9[cps$demState %in% c(26, 35, 42, 28, 24, 17, 16)] <- 4
cps$div9[cps$demState %in% c(8, 9, 21, 47, 49, 41, 34, 10, 11)] <- 5
cps$div9[cps$demState %in% c(18, 43, 25, 1)] <- 6
cps$div9[cps$demState %in% c(4, 19, 37, 44)] <- 7
cps$div9[cps$demState %in% c(13, 27, 51, 29, 45, 6, 3, 32)] <- 8
cps$div9[cps$demState %in% c(2, 48, 38, 12, 5)] <- 9
levels(cps$div9) <- c("New England", "Mid-Atlantic", "East North Central", "West North Central",
                       "South Atlantic", "East South Central", "West South Central", "Mountain", 
                       "Pacific")
table(cps$div9)




################
#-----------Create Weighted Proportion tables for target parameters
################

library(foreign)
library(car)
library(questionr)
require(survey)

#National Sample
#Define target parameters
demos <- c("demRaceFull",   "demGender", "ageDem", "edu5", "demInc", "div9")
#Loop through parameters to create weighted tables
for (i in seq_along(demos)) {    #for each population parameter (demos)
  x <- cps[, demos[i]]    #identify the column with name cps$[demos/variable name]
  a <- wtd.table(x, weights = cps$wts_cps)  #construct a weighted table of counts
  a <- a / sum(a) #Turn counts into proportions
  a <- data.frame(a) #turn list into dataframe
    a0 <- cbind(rep(demos[i], nlevels(as.factor(x))), a)  #identify each newly constructed table with the correspondeding variable it comes from / demos
  if (i == 1)
    b <- a0
  else
    b <- rbind(b, a0)
}
b
#Clean Table
colnames(b) <- c("Item",  "Value", "Pop_Prop")  #rename columns
#b$Source <- "2012CPS"
#class(b$Population_Proportion)
#b$Population_Proportion <- as.numeric(b$Population_Proportion)
cps_nat_parameters <- b
cps_nat_parameters$Pop_Count <- cps_nat_parameters$Pop_Prop*  nrow(data4)
cps_nat_parameters <- cps_nat_parameters[-3]
cps_nat_parameters


#Creat seperate tables
cps_nat_parameters
x <- split(cps_nat_parameters, cps_nat_parameters$Item)
x

race.dist <- x$demRaceFull
race.dist <- race.dist[-1]
colnames(race.dist) <- c("demRaceFull", "Freq")
race.dist

gen.dist <- x$demGender
gen.dist <- gen.dist[-1]
colnames(gen.dist) <- c("demGender", "Freq")
gen.dist

age.dist <- x$ageDem
age.dist <- age.dist[-1]
colnames(age.dist) <- c("ageDem", "Freq")
age.dist

edu.dist <- x$edu5
edu.dist <- edu.dist[-1]
colnames(edu.dist) <- c("edu5", "Freq")
edu.dist

inc.dist <- x$demInc
inc.dist <- inc.dist[-1]
colnames(inc.dist) <- c("demInc", "Freq")
inc.dist

div.dist <- x$div9
div.dist <- div.dist[-1]
colnames(div.dist) <- c("div9", "Freq")
div.dist

#Construct national sample weight
nat.unweight <- svydesign(ids=~1, data=data4)
nat.rake  <- rake(design = nat.unweight,
                   sample.margins = list(~demRaceFull,  ~demGender, ~ageDem, ~edu5, ~demInc, ~div9), 
                   population.margins = list(race.dist, gen.dist, age.dist, edu.dist, inc.dist, div.dist))
summary(weights(nat.rake))

#trim weight
nat <- trimWeights(nat.rake, 
				lower= 0.5,
				upper=6,
                  strict=TRUE)
summary(weights(nat))


#Attach weight to the sample
data4$wts_nat <- weights(nat)
head(data4$wts_nat)







################
#-----------ANALYSIS - MAKE TABLES AND FIGURES FOR STEP 2
################


#Slim down file to just Black, Hispanc, White

# Subset Sample to only Black, White Latino
data5 <- subset(data4, demRaceFull %in% c(1, 2, 6))
table(data5$demRaceFull)

data <- data5
names(data)

# Change order and labels of race3 factors for anlalysis
table(data$demRaceFull)
data$race4 <- factor(data$demRaceFull, levels = c(1, 2, 6), labels = c("White", "Black", "Latino") )
table(data$race4)




#----Calculate weighted mean and 95% confidence intervals
#boot function for calculating weighted mean
sample.wtd.mean <- function(x, wb, d) {
    return(weighted.mean(x = x[d], w = wb[d], na.rm=T ))
} 


#Vote
Like.Vote <- ddply(data, "race4", summarise,
	        mean = weighted.mean(x=vlike, w=wts_nat, na.rm=T),
	        se = sd(boot(vlike, sample.wtd.mean, R = 1000, wb = wts_nat)$t),		
			ci95_t = quantile( boot(vlike, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.025)),
			ci95_b = quantile( boot(vlike, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.975)),		ci80_t = quantile( boot(vlike, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.1), na.rm=T),
		ci80_b = quantile( boot(vlike, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.9), na.rm=T)	
	)		
Like.Vote
#Remove final row of dataframe
Like.vote <-Like.Vote[-4,]
Like.vote


#Rally
Like.Rally <- ddply(.data = data, "race4", .fun = summarise,
        mean = weighted.mean(x=rlike, w=wts_nat, na.rm=T),
        se = sd(boot(rlike, sample.wtd.mean, R = 1000, wb = wts_nat)$t),		
		ci95_t = quantile( boot(rlike, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.025)),
		ci95_b = quantile( boot(rlike, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.975)),		ci80_t = quantile( boot(rlike, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.1), na.rm=T),
		ci80_b = quantile( boot(rlike, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.9), na.rm=T)		
	)		
Like.Rally
#Remove final row of dataframe
Like.rally<- Like.Rally[-4,]
Like.rally


#Potluck
Like.Pot <- ddply(.data = data, .variables= .(race4), .fun = summarise,
        mean = weighted.mean(x=blike, w=wts_nat, na.rm=T),
        se = sd(boot(blike, sample.wtd.mean, R = 1000, wb = wts_nat)$t),		
		ci95_t = quantile( boot(blike, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.025)),
		ci95_b = quantile( boot(blike, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.975)),		ci80_t = quantile( boot(blike, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.1), na.rm=T),
		ci80_b = quantile( boot(blike, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.9), na.rm=T)	
	)		
Like.Pot
#Remove final row of dataframe
Like.pot<- Like.Pot[-4,]
Like.pot





# ------ APPENDIX FIGURE A1

#-------- Plotting Raw Data
#Voter Likability and Potluck
Like.pot$Treat <- "Potluck"
Like.vote$Treat <- "Voter"
Like.rally$Treat <- "Rally"
like.combo.vote <- rbind(Like.pot, Like.vote)
like.combo.vote

# ---- Appendix Figure 1A (a)
FigA1a <- ggplot(like.combo.vote, aes(x = Treat, y = mean)) + 
   geom_point(aes(colour = race4)) + 
			# geom_segment( aes(y = ci95_b, yend = ci95_t,  x = factor(Treat), xend=factor(Treat)) +
 geom_line(aes(colour= race4, group=race4, linetype= race4)) + 
 #scale_colour_manual(values = c( "chartreuse4", "midnightblue" )) + #Set colors of points
  	scale_color_brewer(palette = "Set1")+	
  	 scale_linetype_manual(values =c(2, 1, 3)) +
        theme_bw() + #theme(axis.title.x = theme_text(size = 12, vjust = .25))+ 
#        xlab("Mean Standardized by Potluck Score") +  ylab("Respondent's Race")  +
        xlab(" ") +  ylab("Mean Likability Score ")  +
	scale_y_continuous(limit = c(4.8, 5.6)) +
	  theme(legend.position = "top",   #Move legend to top
   legend.title=element_blank(), 	# Remove legend name
   legend.key=element_blank(), 	#Remove box around legend points
      legend.text = element_text(size = 12), #increase size of legnd text
   legend.key.width=unit(4,"line"), #increase space between legend symbols
   axis.title.x = element_text(vjust = -0.3 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5),	#move x label away from plot
              axis.text=element_text(size=14) #Increase size of axis text
) 
FigA1a

#Save
pdf('Raw_Like_Voter_Labels.pdf')
FigA1a
dev.off()


# ---- Appendix Figure 1A (c)
#------RALLY: Likability Raw Plots
like.combo.ral <- rbind(Like.pot, Like.rally)
like.combo.ral

#Plot: For Paper (Labels, square)
FigA1c <- ggplot(like.combo.ral, aes(x = Treat, y = mean)) + 
   geom_point(aes(colour = race4)) + 
			# geom_segment( aes(y = ci95_b, yend = ci95_t,  x = factor(Treat), xend=factor(Treat)) +
 geom_line(aes(colour= race4, group=race4, linetype= race4)) + 
 #scale_colour_manual(values = c( "chartreuse4", "midnightblue" )) + #Set colors of points
  	scale_color_brewer(palette = "Set1")+	
  	 scale_linetype_manual(values =c(2, 1, 3)) +
        theme_bw() + #theme(axis.title.x = theme_text(size = 12, vjust = .25))+ 
#        xlab("Mean Standardized by Potluck Score") +  ylab("Respondent's Race")  +
        xlab(" ") +  ylab("Mean Likability Score ")  +
	scale_y_continuous(limit = c(4.8, 5.6)) +
	  theme(legend.position = "top",   #Move legend to top
   legend.title=element_blank(), 	# Remove legend name
   legend.key=element_blank(), 	#Remove box around legend points
      legend.text = element_text(size = 12), #increase size of legnd text
   legend.key.width=unit(4,"line"), #increase space between legend symbols
   axis.title.x = element_text(vjust = -0.3 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5),	#move x label away from plot
              axis.text=element_text(size=14) #Increase size of axis text
) 
FigA1c

#Save
pdf('Raw_Like_Rally_Labels.pdf')
FigA1c
dev.off()



# -- Respectability Estimates
#Vote
Resp.Vote <- ddply(.data = data, .variables= .(race4), .fun = summarise,
        mean = weighted.mean(x=v.resp, w=wts_nat, na.rm=T),
        se = sd(boot(v.resp, sample.wtd.mean, R = 1000, wb = wts_nat)$t),
		ci95_t = quantile( boot(v.resp, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.025)),
		ci95_b = quantile( boot(v.resp, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.975)),
			ci80_t = quantile( boot(v.resp, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.1), na.rm=T),
		ci80_b = quantile( boot(v.resp, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.9), na.rm=T)	
		)
Resp.Vote
#Remove final row of dataframe
Resp.vote <-Resp.Vote[-4,]
Resp.vote

#Rally
Resp.Rally <- ddply(.data = data, .variables= .(race4), .fun = summarise,
        mean = weighted.mean(x=r.resp, w=wts_nat, na.rm=T),
        se = sd(boot(r.resp, sample.wtd.mean, R = 1000, wb = wts_nat)$t),
		ci95_t = quantile( boot(r.resp, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.025)),
		ci95_b = quantile( boot(r.resp, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.975)),
			ci80_t = quantile( boot(r.resp, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.1), na.rm=T),
		ci80_b = quantile( boot(r.resp, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.9), na.rm=T)	
		)
Resp.Rally
#Remove final row of dataframe
Resp.rally <-Resp.Rally[-4,]
Resp.rally


#Potluck
Resp.Pot <- ddply(.data = data, .variables= .(race4), .fun = summarise,
        mean = weighted.mean(x=b.resp, w=wts_nat, na.rm=T),
        se = sd(boot(b.resp, sample.wtd.mean, R = 1000, wb = wts_nat)$t),
		ci95_t = quantile( boot(b.resp, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.025)),
		ci95_b = quantile( boot(b.resp, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.975)),
			ci80_t = quantile( boot(b.resp, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.1), na.rm=T),
		ci80_b = quantile( boot(b.resp, sample.wtd.mean, R=1000, wb= wts_nat)$t, c(0.9), na.rm=T)	
		)
Resp.Pot
#Remove final row of dataframe
Resp.pot <-Resp.Pot[-4,]
Resp.pot



#-------- Plotting Raw Data: Voter Respectability
#Voter Likability and Potluck
Resp.pot$Treat <- "Potluck"
Resp.vote$Treat <- "Voter"
resp.combo.vote <- rbind(Resp.pot, Resp.vote)
resp.combo.vote



# ---- Appendix Figure A1b

FigA1b <- ggplot(resp.combo.vote, aes(x = Treat, y = mean)) + 
   geom_point(aes(colour = race4)) + 
			# geom_segment( aes(y = ci95_b, yend = ci95_t,  x = factor(Treat), xend=factor(Treat)) +
 geom_line(aes(colour= race4, group=race4, linetype= race4)) + 
 #scale_colour_manual(values = c( "chartreuse4", "midnightblue" )) + #Set colors of points
  	scale_color_brewer(palette = "Set1")+	
  	 scale_linetype_manual(values =c(2, 1, 3)) +
        theme_bw() + #theme(axis.title.x = theme_text(size = 12, vjust = .25))+ 
#        xlab("Mean Standardized by Potluck Score") +  ylab("Respondent's Race")  +
        xlab(" ") +  ylab("Mean Respectability Score ")  +
	scale_y_continuous(limit = c(4.8, 5.6)) +
	  theme(legend.position = "top",   #Move legend to top
   legend.title=element_blank(), 	# Remove legend name
   legend.key=element_blank(), 	#Remove box around legend points
      legend.text = element_text(size = 12), #increase size of legnd text
   legend.key.width=unit(4,"line"), #increase space between legend symbols
   axis.title.x = element_text(vjust = -0.3 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5),	#move x label away from plot
              axis.text=element_text(size=14) #Increase size of axis text
) 
FigA1b

#Save
pdf('Raw_Resp_Voter_Labels.pdf')
FigA1b
dev.off()


#-------- Plotting Raw Data: Rally Respectability
#Voter Likability and Potluck
Resp.pot$Treat <- "Potluck"
Resp.rally$Treat <- "Rally"
resp.combo.ral <- rbind(Resp.pot, Resp.rally)
resp.combo.ral


# ---- Appendix Figure A1d

FigA1d <- ggplot(resp.combo.ral, aes(x = Treat, y = mean)) + 
   geom_point(aes(colour = race4)) + 
			# geom_segment( aes(y = ci95_b, yend = ci95_t,  x = factor(Treat), xend=factor(Treat)) +
 geom_line(aes(colour= race4, group=race4, linetype= race4)) + 
 #scale_colour_manual(values = c( "chartreuse4", "midnightblue" )) + #Set colors of points
  	scale_color_brewer(palette = "Set1")+	
  	 scale_linetype_manual(values =c(2, 1, 3)) +
        theme_bw() + #theme(axis.title.x = theme_text(size = 12, vjust = .25))+ 
#        xlab("Mean Standardized by Potluck Score") +  ylab("Respondent's Race")  +
        xlab(" ") +  ylab("Mean Respectability Score ")  +
	scale_y_continuous(limit = c(4.8, 5.6)) +
	  theme(legend.position = "top",   #Move legend to top
   legend.title=element_blank(), 	# Remove legend name
   legend.key=element_blank(), 	#Remove box around legend points
      legend.text = element_text(size = 12), #increase size of legnd text
   legend.key.width=unit(4,"line"), #increase space between legend symbols
   axis.title.x = element_text(vjust = -0.3 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5),	#move x label away from plot
              axis.text=element_text(size=14) #Increase size of axis text
) 
FigA1d

#Save
pdf('Raw_Resp_Rally_Labels.pdf')
FigA1d
dev.off()







########
# ----- Make Figures 1 & 2
#######

#Make Vector of Standardized Values - treatment minus control
data$vstan <- data$vlike - data$blike
data$rstan <- data$rlike - data$blike
data$vstan.r <- data$v.resp - data$b.resp
data$rstan.r <- data$r.resp - data$b.resp



#------- DIFFERENCE IN DIFFERENCE ESTIMATES
#---------Compared to Whites, With Plots

#VOTER

#Voter Likability
likevote.dd <- lm(vstan ~ race4, data = data, weight = wts_nat)  #Create difference in difference means
summary(likevote.dd)
c2 <- confint(likevote.dd) #Confidence intervals from regression
likevote.dd2 <- as.data.frame(coef(likevote.dd)) #put coefficients output into a data frame
likevote.dd2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
likevote.dd2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
likevote.dd2$measure <- "Likability" #make a measure indicator
likevote.dd2$cat <-"Aggregate"
likevote.dd2$race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race indicator
likevote.dd2 <- rename(likevote.dd2, c("coef(likevote.dd)" = "mean")) #rename column
likevote.dd2 <- likevote.dd2[-1,] #remove white intercept
likevote.dd2

#Voter Respectability
respvote.dd <- lm(vstan.r ~ race4, data = data, weight = wts_nat)
summary(respvote.dd)
c1<- confint(respvote.dd)
respvote.dd2 <- as.data.frame(coef(respvote.dd)) #put coefficients output into a data frame
respvote.dd2 $ci95_b <- c1[,1]# Put confidence intervals into the data frame 
respvote.dd2 $ci95_t <- c1[,2]# Put confidence intervals into the data frame 
respvote.dd2 $measure <- "Respectability" #make a measure indicator
respvote.dd2$cat <-"Aggregate"
respvote.dd2 $race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race indicator
respvote.dd2 <- rename(respvote.dd2, c("coef(respvote.dd)" = "mean")) 
respvote.dd2 <- respvote.dd2[-1,] #remove white intercept
respvote.dd2


#Create a new low and high ses categories
#Binary SES code
data$ses_bi <- rep(NA, length(data$caseid))
data$ses_bi[data$inc_sm ==1 & data$edu_sm ==1]<-1
data$ses_bi[data$inc_sm ==2 & data$edu_sm ==2]<-2
data$ses_bi[data$inc_sm ==3 & data$edu_sm ==3]<-2
table(data$ses_bi)
table(data$ses_bi, data$race4)
#Low class (1): Income under 30k and High school graduatuion or lesst
#Middle / Upper Class (2): Income ober 30k and At least some college


#SES: Binary
#Voter Likability
voter.like.ses_bi1 <- lm(data$vstan[data$ses_bi ==1]~ data$race4[data$ses_bi ==1], weight = data$wts_nat[data$ses_bi ==1])
summary(voter.like.ses_bi1)
c2 <- confint(voter.like.ses_bi1) #Confidence intervals from regression
voter.like.ses_bi1.2 <- as.data.frame(coef(voter.like.ses_bi1)) #put coefficients output into a data frame
voter.like.ses_bi1.2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
voter.like.ses_bi1.2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
voter.like.ses_bi1.2 $measure <- "Likability" #make a measure indicator
voter.like.ses_bi1.2$cat <- "Low SES"
voter.like.ses_bi1.2 $race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
voter.like.ses_bi1.2 <- rename(voter.like.ses_bi1.2, c("coef(voter.like.ses_bi1)" = "mean")) #rename column
voter.like.ses_bi1.2 <- voter.like.ses_bi1.2[-1,] #remove white intercept
voter.like.ses_bi1.2

#Voter Respectability
voter.resp.ses_bi1 <- lm(data$vstan.r[data$ses_bi ==1]~ data$race4[data$ses_bi ==1], weight = data$wts_nat[data$ses_bi ==1])
summary(voter.resp.ses_bi1)
c2 <- confint(voter.resp.ses_bi1) #Confidence intervals from regression
voter.resp.ses_bi1.2 <- as.data.frame(coef(voter.resp.ses_bi1)) #put coefficients output into a data frame
voter.resp.ses_bi1.2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
voter.resp.ses_bi1.2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
voter.resp.ses_bi1.2 $measure <- "Respectability" #make a measure indicator
voter.resp.ses_bi1.2$cat <- "Low SES"
voter.resp.ses_bi1.2 $race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
voter.resp.ses_bi1.2 <- rename(voter.resp.ses_bi1.2, c("coef(voter.resp.ses_bi1)" = "mean")) #rename column
voter.resp.ses_bi1.2 <- voter.resp.ses_bi1.2[-1,] #remove white intercept
voter.resp.ses_bi1.2 <- voter.resp.ses_bi1.2[-c(3:8),] #remove white intercept
voter.resp.ses_bi1.2

voter.ses_bi1.dd <- rbind(voter.like.ses_bi1.2, voter.resp.ses_bi1.2)
voter.ses_bi1.dd

#AMONG SES = High/Middle (2)
#Voter Likability
voter.like.ses_bi2 <- lm(data$vstan[data$ses_bi ==2]~ data$race4[data$ses_bi ==2], weight = data$wts_nat[data$ses_bi ==2])
summary(voter.like.ses_bi2)
c2 <- confint(voter.like.ses_bi2) #Confidence intervals from regression
voter.like.ses_bi2.2 <- as.data.frame(coef(voter.like.ses_bi2)) #put coefficients output into a data frame
voter.like.ses_bi2.2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
voter.like.ses_bi2.2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
voter.like.ses_bi2.2 $measure <- "Likability" #make a measure indicator
voter.like.ses_bi2.2$cat <- "Med/High SES"
voter.like.ses_bi2.2 $race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
voter.like.ses_bi2.2 <- rename(voter.like.ses_bi2.2, c("coef(voter.like.ses_bi2)" = "mean")) #rename column
voter.like.ses_bi2.2 <- voter.like.ses_bi2.2[-1,] #remove white intercept
voter.like.ses_bi2.2

#Voter Respectability
voter.resp.ses_bi2 <- lm(data$vstan.r[data$ses_bi ==2]~ data$race4[data$ses_bi ==2], weight = data$wts_nat[data$ses_bi ==2])
summary(voter.resp.ses_bi2)
c2 <- confint(voter.resp.ses_bi2) #Confidence intervals from regression
voter.resp.ses_bi2.2 <- as.data.frame(coef(voter.resp.ses_bi2)) #put coefficients output into a data frame
voter.resp.ses_bi2.2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
voter.resp.ses_bi2.2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
voter.resp.ses_bi2.2 $measure <- "Respectability" #make a measure indicator
voter.resp.ses_bi2.2$cat <- "Med/High SES"
voter.resp.ses_bi2.2 $race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
voter.resp.ses_bi2.2 <- rename(voter.resp.ses_bi2.2, c("coef(voter.resp.ses_bi2)" = "mean")) #rename column
voter.resp.ses_bi2.2 <- voter.resp.ses_bi2.2[-1,] #remove white intercept
voter.resp.ses_bi2.2 <- voter.resp.ses_bi2.2[-c(3:8),] #remove white intercept
voter.resp.ses_bi2.2


#ALONG GENDER
#Among Women
table(data$gender)
#Voter Likability
#Work with subsetting gender data
fem <- subset(data, gender=="Female")
table(fem$gender)
voter.like.w <- lm(vstan~ race4, data = fem, weight = wts_nat)
summary(voter.like.w)
c2 <- confint(voter.like.w) #Confidence intervals from regression
voter.like.w2 <- as.data.frame(coef(voter.like.w)) #put coefficients output into a data frame
voter.like.w2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
voter.like.w2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
voter.like.w2$measure <- "Likability" #make a measure indicator
voter.like.w2$cat <- "Women"
voter.like.w2$race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
voter.like.w2 <- rename(voter.like.w2, c("coef(voter.like.w)" = "mean")) #rename column
voter.like.w2 <- voter.like.w2[-1,] #remove white intercept
voter.like.w2 <- voter.like.w2[-c(3:8),] #remove white intercept
voter.like.w2

#Voter Respectability
voter.res.w <- lm(vstan.r~ race4, data = fem, weight = wts_nat)
summary(voter.res.w)
c2 <- confint(voter.res.w) #Confidence intervals from regression
voter.res.w2 <- as.data.frame(coef(voter.res.w)) #put coefficients output into a data frame
voter.res.w2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
voter.res.w2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
voter.res.w2 $measure <- "Respectability" #make a measure indicator
voter.res.w2$cat <- "Women"
voter.res.w2 $race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
voter.res.w2 <- rename(voter.res.w2, c("coef(voter.res.w)" = "mean")) #rename column
voter.res.w2 <- voter.res.w2[-1,] #remove white intercept
voter.res.w2 <- voter.res.w2[-c(3:8),] #remove white intercept
voter.res.w2


#Among Men
table(data$gender)
#Voter Likability
#Work with subsetting gender data
mal <- subset(data, gender=="Male")
table(mal$gender)
voter.like.m <- lm(vstan~ race4, data = mal, weight = wts_nat)
summary(voter.like.m)
c2 <- confint(voter.like.m) #Confidence intervals from regression
voter.like.m2 <- as.data.frame(coef(voter.like.m)) #put coefficients output into a data frame
voter.like.m2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
voter.like.m2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
voter.like.m2$measure <- "Likability" #make a measure indicator
voter.like.m2$cat <- "Men"
voter.like.m2$race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
voter.like.m2 <- rename(voter.like.m2, c("coef(voter.like.m)" = "mean")) #rename column
voter.like.m2 <- voter.like.m2[-1,] #remove mhite intercept
voter.like.m2 <- voter.like.m2[-c(3:8),] #remove mhite intercept
voter.like.m2

#Voter Respectability
voter.res.m <- lm(vstan.r~ race4, data = mal, weight = wts_nat)
summary(voter.res.m)
c2 <- confint(voter.res.m) #Confidence intervals from regression
voter.res.m2 <- as.data.frame(coef(voter.res.m)) #put coefficients output into a data frame
voter.res.m2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
voter.res.m2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
voter.res.m2 $measure <- "Respectability" #make a measure indicator
voter.res.m2$cat <- "Men"
voter.res.m2 $race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
voter.res.m2 <- rename(voter.res.m2, c("coef(voter.res.m)" = "mean")) #rename column
voter.res.m2 <- voter.res.m2[-1,] #remove white intercept
voter.res.m2 <- voter.res.m2[-c(3:8),] #remove white intercept
voter.res.m2


# Bind together Estimates in a single dataset
Joined <- rbind(likevote.dd2, respvote.dd2, voter.like.ses_bi1.2, voter.resp.ses_bi1.2, voter.like.ses_bi2.2, voter.resp.ses_bi2.2, voter.like.w2, voter.res.w2, voter.like.m2, voter.res.m2)
names(Joined)
Joined


#Subset to Black respondents
join_black <- subset(Joined, race4=="Black Respondents")
join_black
#Aggregate Black 
join_black_agg_v <- join_black [c(1,2),]
join_black_agg_v
#Remove aggregate: subplots 
join_black_sub_v <- join_black[-c(1,2),]

#Subset to Lation respondents
join_latino <- subset(Joined, race4=="Latino Respondents")
join_latino
#Aggregate Latino 
join_latino_agg_v <- join_latino [c(1,2),]
join_latino_agg_v
#Subsetted Latino Groups
join_latino_sub_v<- join_latino[-c(1,2),]

#Join Black and Latino Aggregage Results
join_agg_v <- rbind(join_black_agg_v, join_latino_agg_v)
join_agg_v

#Reorder levels for plotting
#Black
join_black_sub_v$cat <- factor(join_black_sub_v$cat, levels = c("Low SES", "Med/High SES", "Women", "Men"))
join_black_sub_v
#Latino
join_latino_sub_v$cat <- factor(join_latino_sub_v$cat, levels = c( "Low SES", "Med/High SES", "Women", "Men"))
join_latino_sub_v



#------ Figure 1 (b)
#Plot -- Aggregate Results 
Fig1b <- ggplot(join_agg_v, aes(x = factor(measure), y = mean)) + 
        geom_point( ) + 
		geom_errorbar( aes(ymin = ci95_b, ymax= ci95_t), width = 0.1) +
        facet_wrap(~race4, ncol=5) +
#scale_colour_manual(values = c( "#377eb8", "#4daf4a" )) + #Set colors of points

        theme_bw() +   #theme(axis.title.x = theme_text(size = 12, vjust = .25))+ 
#        xlab("Mean Standardized by Potluck Score") +  ylab("Respondent's Race")  +
        xlab(" ") +  ylab("Estimated Effect of Voting")  +
		scale_y_continuous(limit = c(-0.2, 0.61), breaks=c(-0.4, -0.2, 0, 0.2,  0.4, 0.6 )) + 
		scale_x_discrete(breaks=c("Likability", "Respectability"),labels=c("Like", "Respect")) +
	geom_hline(yintercept=0, slope=0, linetype=2, colour = "#e41a1c") +
	theme(legend.position="none",
		strip.text = element_text(size=20),
	             axis.text=element_text(size=17), #Increase size of axis text
   axis.title.x = element_text(vjust = -0.3, size = 17 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5, size = 17)	#move x label away from plot
	)
Fig1b
#Save
pdf( "Aggregate_Voting.pdf", width=8, height=8)
Fig1b
dev.off()



#Figure 2(a)
#Plot -- Black Across Groups
Fig2a <- ggplot(join_black_sub_v, aes(x = factor(measure), y = mean)) + 
        geom_point( ) + 
		geom_errorbar( aes(ymin = ci95_b, ymax= ci95_t), width = 0.1) +
        facet_wrap(~cat, ncol=5) +
#scale_colour_manual(values = c( "#377eb8", "#4daf4a" )) + #Set colors of points

        theme_bw() +   #theme(axis.title.x = theme_text(size = 12, vjust = .25))+ 
#        xlab("Mean Standardized by Potluck Score") +  ylab("Respondent's Race")  +
        xlab(" ") +  ylab("Estimated Effect of Voting")  +
		scale_y_continuous(limit = c(-0.4, 0.81), breaks=c(-0.4, -0.2, 0, 0.2,  0.4, 0.6, 0.8 )) + 
		scale_x_discrete(breaks=c("Likability", "Respectability"),labels=c("Like", "Respect")) +
	geom_hline(yintercept=0, slope=0, linetype=2, colour = "#e41a1c") +
	theme(legend.position="none",
		strip.text = element_text(size=20),
	             axis.text=element_text(size=17), #Increase size of axis text
   axis.title.x = element_text(vjust = -0.3, size = 17 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5, size = 17)	#move x label away from plot
	)
Fig2a

#Save
pdf( "Black_Subset.pdf", width=15, height=4)
Fig2a
dev.off()


# Figure 2(c)
#Plot -- Latino Across Groups
Fig2c <- ggplot(join_latino_sub_v, aes(x = factor(measure), y = mean)) + 
        geom_point( ) + 
		geom_errorbar( aes(ymin = ci95_b, ymax= ci95_t), width = 0.1) +
        facet_wrap(~cat, ncol=5) +
        theme_bw() +   #theme(axis.title.x = theme_text(size = 12, vjust = .25))+ 
        xlab(" ") +  ylab("Estimated Effect of Voting")  +
		scale_y_continuous(limit = c(-0.4, 0.81), breaks=c(-0.4, -0.2, 0, 0.2,  0.4, 0.6, 0.8 )) + 
		scale_x_discrete(breaks=c("Likability", "Respectability"),labels=c("Like", "Respect")) +
	geom_hline(yintercept=0, slope=0, linetype=2, colour = "#e41a1c") +
	theme(legend.position="none",
		strip.text = element_text(size=20),
	             axis.text=element_text(size=17), #Increase size of axis text
   axis.title.x = element_text(vjust = -0.3, size = 17 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5, size = 17)	#move x label away from plot
	)
Fig2c

#Save
pdf( "Latino_Subset.pdf", width=15, height=4)
Fig2c
dev.off()


#------- RALLY DIF IN DIF PLOTS
#rallyr Likability
likerally.dd <- lm(rstan ~ race4, data = data, weight = wts_nat)  #Create difference in difference means
summary(likerally.dd)
c2 <- confint(likerally.dd) #Confidence intervals from regression
likerally.dd2 <- as.data.frame(coef(likerally.dd)) #put coefficients output into a data frame
likerally.dd2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
likerally.dd2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
likerally.dd2$measure <- "Likability" #make a measure indicator
likerally.dd2$cat <-"Aggregate"
likerally.dd2$race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race indicator
likerally.dd2 <- rename(likerally.dd2, c("coef(likerally.dd)" = "mean")) #rename column
likerally.dd2 <- likerally.dd2[-1,] #remove white intercept
likerally.dd2

#rallyr Respectability
resprally.dd <- lm(rstan.r ~ race4, data = data, weight = wts_nat)
summary(resprally.dd)
c1<- confint(resprally.dd)
resprally.dd2 <- as.data.frame(coef(resprally.dd)) #put coefficients output into a data frame
resprally.dd2 $ci95_b <- c1[,1]# Put confidence intervals into the data frame 
resprally.dd2 $ci95_t <- c1[,2]# Put confidence intervals into the data frame 
resprally.dd2 $measure <- "Respectability" #make a measure indicator
resprally.dd2$cat <-"Aggregate"
resprally.dd2 $race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race indicator
resprally.dd2 <- rename(resprally.dd2, c("coef(resprally.dd)" = "mean")) 
resprally.dd2 <- resprally.dd2[-1,] #remove white intercept
resprally.dd2


#Women
#rally Likability
rallyr.like.wom <- lm(rstan~ race4, data = fem, weight = wts_nat)
summary(rallyr.like.wom)
c2 <- confint(rallyr.like.wom) #Confidence intervals from regression
rallyr.like.wom2 <- as.data.frame(coef(rallyr.like.wom)) #put coefficients output into a data frame
rallyr.like.wom2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
rallyr.like.wom2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
rallyr.like.wom2$measure <- "Likability" #make a measure indicator
rallyr.like.wom2$cat <- "Women"
rallyr.like.wom2$race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
rallyr.like.wom2 <- rename(rallyr.like.wom2, c("coef(rallyr.like.wom)" = "mean")) #rename column
rallyr.like.wom2 <- rallyr.like.wom2[-1,] #remove white intercept
rallyr.like.wom2 <- rallyr.like.wom2[-c(3:8),] #remove white intercept
rallyr.like.wom2

#rallyr Respectability
rallyr.res.wom <- lm(rstan.r~ race4, data = fem, weight = wts_nat)
summary(rallyr.res.wom)
c2 <- confint(rallyr.res.wom) #Confidence intervals from regression
rallyr.res.wom2 <- as.data.frame(coef(rallyr.res.wom)) #put coefficients output into a data frame
rallyr.res.wom2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
rallyr.res.wom2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
rallyr.res.wom2 $measure <- "Respectability" #make a measure indicator
rallyr.res.wom2$cat <- "Women"
rallyr.res.wom2 $race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
rallyr.res.wom2 <- rename(rallyr.res.wom2, c("coef(rallyr.res.wom)" = "mean")) #rename column
rallyr.res.wom2 <- rallyr.res.wom2[-1,] #remove white intercept
rallyr.res.wom2 <- rallyr.res.wom2[-c(3:8),] #remove white intercept
rallyr.res.wom2


#Men
#rally Likability
rallyr.like.men <- lm(rstan~ race4, data = mal, weight = wts_nat)
summary(rallyr.like.men)
c2 <- confint(rallyr.like.men) #Confidence intervals from regression
rallyr.like.men2 <- as.data.frame(coef(rallyr.like.men)) #put coefficients output into a data frame
rallyr.like.men2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
rallyr.like.men2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
rallyr.like.men2$measure <- "Likability" #make a measure indicator
rallyr.like.men2$cat <- "Men"
rallyr.like.men2$race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
rallyr.like.men2 <- rename(rallyr.like.men2, c("coef(rallyr.like.men)" = "mean")) #rename column
rallyr.like.men2 <- rallyr.like.men2[-1,] #remove white intercept
rallyr.like.men2 <- rallyr.like.men2[-c(3:8),] #remove white intercept
rallyr.like.men2

#rallyr Respectability
rallyr.res.men <- lm(rstan.r~ race4, data = mal, weight = wts_nat)
summary(rallyr.res.men)
c2 <- confint(rallyr.res.men) #Confidence intervals from regression
rallyr.res.men2 <- as.data.frame(coef(rallyr.res.men)) #put coefficients output into a data frame
rallyr.res.men2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
rallyr.res.men2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
rallyr.res.men2 $measure <- "Respectability" #make a measure indicator
rallyr.res.men2$cat <- "Men"
rallyr.res.men2 $race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
rallyr.res.men2 <- rename(rallyr.res.men2, c("coef(rallyr.res.men)" = "mean")) #rename column
rallyr.res.men2 <- rallyr.res.men2[-1,] #remove white intercept
rallyr.res.men2 <- rallyr.res.men2[-c(3:8),] #remove white intercept
rallyr.res.men2



#SES Binary
#AMONG SES = 1
#rallyr Likability
rallyr.like.ses_bi1 <- lm(data$rstan[data$ses_bi ==1]~ data$race4[data$ses_bi ==1], weight = data$wts_nat[data$ses_bi ==1])
summary(rallyr.like.ses_bi1)
c2 <- confint(rallyr.like.ses_bi1) #Confidence intervals from regression
rallyr.like.ses_bi1.2 <- as.data.frame(coef(rallyr.like.ses_bi1)) #put coefficients output into a data frame
rallyr.like.ses_bi1.2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
rallyr.like.ses_bi1.2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
rallyr.like.ses_bi1.2 $measure <- "Likability" #make a measure indicator
rallyr.like.ses_bi1.2$cat <- "Low SES"
rallyr.like.ses_bi1.2 $race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
rallyr.like.ses_bi1.2 <- rename(rallyr.like.ses_bi1.2, c("coef(rallyr.like.ses_bi1)" = "mean")) #rename column
rallyr.like.ses_bi1.2 <- rallyr.like.ses_bi1.2[-1,] #remove white intercept
rallyr.like.ses_bi1.2

#rallyr Respectability
rallyr.resp.ses_bi1 <- lm(data$rstan.r[data$ses_bi ==1]~ data$race4[data$ses_bi ==1], weight = data$wts_nat[data$ses_bi ==1])
summary(rallyr.resp.ses_bi1)
c2 <- confint(rallyr.resp.ses_bi1) #Confidence intervals from regression
rallyr.resp.ses_bi1.2 <- as.data.frame(coef(rallyr.resp.ses_bi1)) #put coefficients output into a data frame
rallyr.resp.ses_bi1.2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
rallyr.resp.ses_bi1.2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
rallyr.resp.ses_bi1.2 $measure <- "Respectability" #make a measure indicator
rallyr.resp.ses_bi1.2$cat <- "Low SES"
rallyr.resp.ses_bi1.2 $race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
rallyr.resp.ses_bi1.2 <- rename(rallyr.resp.ses_bi1.2, c("coef(rallyr.resp.ses_bi1)" = "mean")) #rename column
rallyr.resp.ses_bi1.2 <- rallyr.resp.ses_bi1.2[-1,] #remove white intercept
rallyr.resp.ses_bi1.2 <- rallyr.resp.ses_bi1.2[-c(3:8),] #remove white intercept
rallyr.resp.ses_bi1.2

rallyr.ses_bi1.dd <- rbind(rallyr.like.ses_bi1.2, rallyr.resp.ses_bi1.2)
rallyr.ses_bi1.dd

#AMONG SES = 2
#rallyr Likability
rallyr.like.ses_bi2 <- lm(data$rstan[data$ses_bi ==2]~ data$race4[data$ses_bi ==2], weight = data$wts_nat[data$ses_bi ==2])
summary(rallyr.like.ses_bi2)
c2 <- confint(rallyr.like.ses_bi2) #Confidence intervals from regression
rallyr.like.ses_bi2.2 <- as.data.frame(coef(rallyr.like.ses_bi2)) #put coefficients output into a data frame
rallyr.like.ses_bi2.2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
rallyr.like.ses_bi2.2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
rallyr.like.ses_bi2.2 $measure <- "Likability" #make a measure indicator
rallyr.like.ses_bi2.2$cat <- "Med/High SES"
rallyr.like.ses_bi2.2 $race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
rallyr.like.ses_bi2.2 <- rename(rallyr.like.ses_bi2.2, c("coef(rallyr.like.ses_bi2)" = "mean")) #rename column
rallyr.like.ses_bi2.2 <- rallyr.like.ses_bi2.2[-1,] #remove white intercept
rallyr.like.ses_bi2.2

#rallyr Respectability
rallyr.resp.ses_bi2 <- lm(data$rstan.r[data$ses_bi ==2]~ data$race4[data$ses_bi ==2], weight = data$wts_nat[data$ses_bi ==2])
summary(rallyr.resp.ses_bi2)
c2 <- confint(rallyr.resp.ses_bi2) #Confidence intervals from regression
rallyr.resp.ses_bi2.2 <- as.data.frame(coef(rallyr.resp.ses_bi2)) #put coefficients output into a data frame
rallyr.resp.ses_bi2.2 $ci95_b <- c2[,1]# Put confidence intervals into the data frame 
rallyr.resp.ses_bi2.2 $ci95_t <- c2[,2]# Put confidence intervals into the data frame 
rallyr.resp.ses_bi2.2 $measure <- "Respectability" #make a measure indicator
rallyr.resp.ses_bi2.2$cat <- "Med/High SES"
rallyr.resp.ses_bi2.2 $race4 <- c("White Respondents", "Black Respondents", "Latino Respondents") #make a race4 indicator
rallyr.resp.ses_bi2.2 <- rename(rallyr.resp.ses_bi2.2, c("coef(rallyr.resp.ses_bi2)" = "mean")) #rename column
rallyr.resp.ses_bi2.2 <- rallyr.resp.ses_bi2.2[-1,] #remove white intercept
rallyr.resp.ses_bi2.2 <- rallyr.resp.ses_bi2.2[-c(3:8),] #remove white intercept
rallyr.resp.ses_bi2.2



# Bind together Estimates in a single dataset
Joined_rally <- rbind(likerally.dd2, resprally.dd2, 
rallyr.like.ses_bi1.2, rallyr.resp.ses_bi1.2, rallyr.like.ses_bi2.2, rallyr.resp.ses_bi2.2, rallyr.like.wom2, rallyr.res.wom2, rallyr.like.men2, rallyr.res.men2)
names(Joined_rally)
Joined_rally

#Seperate Out Black and Latino Estimates

#Subset to Black respondents
join_black_r <- subset(Joined_rally, race4=="Black Respondents")
join_black_r
#Aggregate Black
join_black_agg <- join_black_r[c(1,2),]
join_black_agg
#Subsetted Black Groups
join_black_sub <- join_black_r[-c(1,2),]

#Subset to Lation respondents
join_latino_r <- subset(Joined_rally, race4=="Latino Respondents")
join_latino_r
#Aggregate Latino 
join_latino_agg <- join_latino_r[c(1,2),]
join_latino_agg
#Subsetted Latino Groups
join_latino_sub<- join_latino_r[-c(1,2),]

#Join Black and Latino Aggregage Results
join_agg <- rbind(join_black_agg, join_latino_agg)
join_agg
#Label Facets for Plotting 
#join_agg$label <- c("Black:White Voting", "Black:White Voting", "Latino:White Voting", "Latino:White Voting") 


# ---- Figure 1A
#Plot -- Black and Latino Aggregate Results
Fig1a <- ggplot(join_agg, aes(x = factor(measure), y = mean)) + 
        geom_point( ) + 
		geom_errorbar( aes(ymin = ci95_b, ymax= ci95_t), width = 0.1) +
        facet_wrap(~race4, ncol=2) +
        theme_bw() +   #theme(axis.title.x = theme_text(size = 12, vjust = .25))+ 
        xlab(" ") +  ylab("Estimated Effect of Rally")  +
		scale_y_continuous(limit = c(-0.2, 0.61), breaks=c(-0.2, 0, 0.2,  0.4, 0.6 )) + 
		scale_x_discrete(breaks=c("Likability", "Respectability"),labels=c("Like", "Respect")) +
	geom_hline(yintercept=0, slope=0, linetype=2, colour = "#e41a1c") +
	theme(legend.position="none",
		strip.text = element_text(size=20),
	             axis.text=element_text(size=17), #Increase size of axis text
   axis.title.x = element_text(vjust = -0.3, size = 17 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5, size = 17)	#move x label away from plot
	)
Fig1a
#Save
pdf("Aggregate_Rally.pdf", width=8, height=8)
Fig1a
dev.off()

#Reorder levels for plotting
#Black
join_black_sub$cat <- factor(join_black_sub$cat, levels = c( "Low SES", "Med/High SES", "Women", "Men"))
join_black_r
#Latino
join_latino_sub$cat <- factor(join_latino_sub$cat, levels = c( "Low SES", "Med/High SES", "Women", "Men"))
join_latino_r



# ----- Figure 2(b)
#Plot -- Black Across Groups
Fig2b <- ggplot(join_black_sub, aes(x = factor(measure), y = mean)) + 
        geom_point( ) + 
		geom_errorbar( aes(ymin = ci95_b, ymax= ci95_t), width = 0.1) +
        facet_wrap(~cat, ncol=5) +
#scale_colour_manual(values = c( "#377eb8", "#4daf4a" )) + #Set colors of points

        theme_bw() +   #theme(axis.title.x = theme_text(size = 12, vjust = .25))+ 
#        xlab("Mean Standardized by Potluck Score") +  ylab("Respondent's Race")  +
        xlab(" ") +  ylab("Estimated Effect of Rally")  +
		scale_y_continuous(limit = c(-0.4, 0.81), breaks=c(-0.4, -0.2, 0, 0.2,  0.4, 0.6, 0.8 )) + 
		scale_x_discrete(breaks=c("Likability", "Respectability"),labels=c("Like", "Respect")) +
	geom_hline(yintercept=0, slope=0, linetype=2, colour = "#e41a1c") +
	theme(legend.position="none",
		strip.text = element_text(size=20),
	             axis.text=element_text(size=17), #Increase size of axis text
   axis.title.x = element_text(vjust = -0.3, size = 17 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5, size = 17)	#move x label away from plot
	)
Fig2b
#Save
pdf( "Black_Subset_Rally.pdf", width=15, height=4)
Fig2b
dev.off()

# --- Figure 2(d)
#Plot -- Latino Across Groups
Fig2d <- ggplot(join_latino_sub, aes(x = factor(measure), y = mean)) + 
        geom_point( ) + 
		geom_errorbar( aes(ymin = ci95_b, ymax= ci95_t), width = 0.1) +
        facet_wrap(~cat, ncol=5) +
        theme_bw() +   #theme(axis.title.x = theme_text(size = 12, vjust = .25))+ 
        xlab(" ") +  ylab("Estimated Effect of Rally")  +
		scale_y_continuous(limit = c(-0.24, 0.85), breaks=c(-0.4, -0.2, 0, 0.2,  0.4, 0.6, 0.8)) + 
		scale_x_discrete(breaks=c("Likability", "Respectability"),labels=c("Like", "Respect")) +
	geom_hline(yintercept=0, slope=0, linetype=2, colour = "#e41a1c") +
	theme(legend.position="none",
		strip.text = element_text(size=20),
	             axis.text=element_text(size=17), #Increase size of axis text
   axis.title.x = element_text(vjust = -0.3, size = 17 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5, size = 17)	#move x label away from plot
	)
Fig2d

#Save
pdf( "Latino_Subset_Rally.pdf", width=15, height=4)
Fig2d
dev.off()





#### --- Data for Appendix Table A2

#  Examining the first appearance versus the total appearance -- checking for period and carryover effects
#Likability
 mean(data$APA002[data$p2act==3]) #Vote first
 mean(data$vlike)  #Vote total
 
 mean(data$APA002[data$p2act==2] ) # Rally first
 mean(data$rlike) #Rally total
 
t.test(data$vlike, data$APA002[data$p2act==3])
t.test(data$rlike, data$APA002[data$p2act==2])

#Respectability
mean(data$APA003[data$p2act==3] ) #Voter first
mean(data$v.resp) #Voter total

mean(data$APA003[data$p2act==2] )
mean(data$r.resp)

t.test(data$APA003[data$p2act==3], data$v.resp)
t.test(data$APA003[data$p2act==2], data$r.resp)







# # #------ Testing Power
# #install.packages("pwr")
# #Info: https://www.statmethods.net/stats/power.html
# require(pwr)

# n<- 1254+1233+744
# bn<- 1233 #calculate more sysnalticallly from data
# wn<- 1254
# ln <- 744

# #For Black-White Test
# pwr.r.test(n=(bn+wn)*2, sig.level=0.05, power=0.8 )
# pwr.r.test(n=(bn+wn)*2, sig.level=0.05, r=0.05)


# #For Latino-White Test
# pwr.r.test(n=(ln+wn)*2, sig.level=0.05, power=0.8 )
# pwr.r.test(n=(ln+wn)*2, sig.level=0.05, r=0.05 )

# #Calculate substantive magnitude on a 1-7 scale
# (0.0397*6)

# (0.19)/(7-1)






############## 
#---------- STEP 3: GEO DATA
#############

data_survey1 <- data
names(data_survey1)

#Add a leading zero to 4 digit zips to match census
library(stringr)
data_survey1$zip <- as.numeric(str_pad(data_survey1$lookupzip, 5, pad = "0"))



#----------Load in Census Zipcode Data: ACS 2009 - 2013
#Nationality and CItizenship data
cen1<- read.csv(file="3ACS_Citizen_Slim.csv")
names(cen1)
length(cen1$ZCTA5A)
#ZCTA5A : zip code 5-digit

#Race, Latino, Language data
cen2<- read.csv(file="3ACS_Race_Slim.csv")
names(cen2)
length(cen2$ZCTA5A)

#Merge together two census files
cen3 <- merge(cen1, cen2, by = "ZCTA5A")
names(cen3)
length(cen3$ZCTA5A)

#Rename zipcode indiactor
cen3$zip <- cen3$ZCTA5A

#merge by zipcode with survey data
data8 <- merge(data_survey1, cen3, by = "zip", all.x = TRUE)
length(data8$caseid)


#Create Proportions for Relevant Categories
#Racial Composition of Zipcode
data8$prop.w <- data8$UEQE002 / data8$UEQE001  #proportion of zipcode that's White 
head(data8$prop.w)
#hist(data8$prop.w)
data8$prop.b <- data8$UEQE003 / data8$UEQE001  #proportion of zipcode that's Black 
head(data8$prop.b)
#hist(data8$prop.b)
data8$prop.h <- data8$UEYE012 / data8$UEYE001  #proportion of zipcode that's hispanic/latino any race
head(data8$prop.h)
#hist(data8$prop.h)

data8$prop.a <- data8$UEQE005/data8$UEQE001 # propoprtion Asian for Robustness checks


#Citizenship
# % of zip code that is US citizens born in the United States
data8$cit <-  data8$UPFE002 / data8$UPFE001
head(data8$cit)


#Create Merged "Average" Respect and Like Measure for DV

data8$v.mean <- (data8$v.resp +data8$vlike)/2
head(data8$v.mean)
head(data8$v.resp)

data8$r.mean <- (data8$r.resp +data8$rlike)/2
head(data8$r.mean)
head(data8$r.resp)

#Reverse Code US Born Citizens
data8$noncit <- 1 - data8$cit


#Normalize to 0-1
#Rally
summary(data8$r.mean)
data8$r.mean.stan <- (data8$r.mean-1)/(7-1)
head(data8$r.mean)
head(data8$r.mean.stan)
#hist(data_LW$r.mean, breaks=6)
#hist(data_LW$r.mean.stan, breaks=7) 
#Voter
summary(data8$v.mean)
data8$v.mean.stan <- (data8$v.mean-1)/(7-1)
head(data8$v.mean)
head(data8$v.mean.stan)


#Re-order Race
data8$race5 <- recode(data8$race4, "'White' = 3; 'Black' = 1; 'Latino' = 2")
table(data8$race5)

#---- Subset Data to Black and White, or Hispanic and White
data_BW <- subset(data8, race4=="Black" | race4=="White")
table(data_BW$race5)
data_BW$race5 <- factor(data_BW$race5, levels=c(1,3), labels=c("Black", "White"))
table(data_BW$race5)

data_LW <- subset(data8, race4=="Latino" | race4=="White")
table(data_LW$race4)
table(data_LW$race5)
data_LW$race5 <- factor(data_LW$race5, levels=c(2,3), labels=c("Latino", "White"))
table(data_LW$race5)



#-------Descriptive Statistics on zipcodes 


summary(data8$UEQE001)
#  Histogram: Showing distribution of zipcode density in total sample
h<- hist(data8$UEQE001, breaks=10, main=" ", xlab="Total Population in Zip Code")
h$density <- h$counts/sum(h$counts)*100


# ---- Figure Appendix A2
pdf("Hist_TotalPopZip_Percent.pdf", width = 8, height = 8)
plot(h, freq=FALSE, ylab="Percent Sample", main=" ", xlab="Total Number of Residents in Zip Code")
 main="Distribution of Population Density by Zip Code in Survey Sample"
dev.off()


require(lattice)

# Appendix Figure A3
pdf("Hist_BlackZip.pdf", width = 8, height = 8)
histogram(~ prop.b | race4, data=data_BW, xlab="Proportion Zip Code Black")
dev.off()

# Appendix Figure A4
pdf("Hist_LatinoZip.pdf", width = 8, height = 8)
histogram(~ prop.h | race4, data=data_LW, xlab="Proportion Zip Code Latino")
dev.off()

# Appendix Figure A5
pdf("Hist_CitZip.pdf", width = 8, height = 8)
histogram(~ noncit | race4, data=data_LW, xlab="Proportion Zip Code Latino")
dev.off()






#----- Regressions and Make Plots
#install.packages("jtools")
require(jtools)


# ----- Make Panels of Figure 3


#-----------Black in Zipcode

data_BW$race4 <- factor(data_BW$race4, levels=c("White","Black"), labels=c("White", "Black"))

#Average Voter -- Black vs White in Increasingly Black Neighborhood
 b1 <- lm(v.mean.stan ~prop.b*race5 + inc + edu1+gender, data=data_BW, weight=wts_nat)
summary(b1)



b1plot <- interact_plot(b1, pred = "prop.b", modx = "race5", interval = TRUE, int.width = 0.8, color.class="Set2", x.label = "Proportion Zip Code Black", y.label="Social Value") + theme_bw() +
  	 #scale_linetype_manual(values =c(2, 1, 3)) +
	scale_y_continuous(breaks=seq(0.56, 0.84, 0.04)) + expand_limits(y=c(0.56,0.84)) +
	  theme(legend.position = "top",   #Move legend to top
   legend.title=element_blank(), 	# Remove legend name
   legend.key=element_blank(), 	#Remove box around legend points
      legend.text = element_text(size = 28), #increase size of legnd text
   legend.key.width=unit(4,"line"), #increase space between legend symbols
             axis.text=element_text(size=25), #Increase size of axis text
   axis.title.x = element_text(vjust = -0.3, size=28 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5, size=28)	#move x label away from plot
) 
b1plot

pdf('Black_Zip_Two_Vote_1plot.pdf', width  = 8, height = 9)
b1plot
dev.off()


 #Average Rally -- Black vs. White in Increasingly Black Neighborhoods
b2 <- lm(r.mean.stan ~prop.b*race5 + inc + edu1+gender, data=data_BW, weight=wts_nat)
summary(b2)
#b2plot<- plot(effect("prop.b:race5", b2),   multiline=FALSE, ylab="Avg. Social Reward for Ralliers", xlab="Proportion Black Residents", rug=TRUE, main=NULL)
b2plot <- interact_plot(b2, pred = "prop.b", modx = "race5", interval = TRUE, int.width = 0.8, color.class="Set2", x.label = "Proportion Zip Code Black", y.label="Social Value") + theme_bw() +
  	 #scale_linetype_manual(values =c(2, 1, 3)) +
	scale_y_continuous(breaks=seq(0.56, 0.84, 0.04)) + expand_limits(y=c(0.56,0.84)) +
	  theme(legend.position = "top",   #Move legend to top
   legend.title=element_blank(), 	# Remove legend name
   legend.key=element_blank(), 	#Remove box around legend points
      legend.text = element_text(size = 28), #increase size of legnd text
   legend.key.width=unit(4,"line"), #increase space between legend symbols
             axis.text=element_text(size=25), #Increase size of axis text
   axis.title.x = element_text(vjust = -0.3, size=28 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5, size=28)	#move x label away from plot
) 
b2plot

pdf('Black_Zip_Two_Rally_1plot.pdf', width  = 8, height = 9)
b2plot
dev.off()


#----------Proportion Citizens
#Average Voter - Latino vs. White by Zip code people either non-us born cits, or non-cits
c1 <- lm(v.mean.stan ~noncit*race5 + inc + edu1 +gender, data=data_LW, weight=wts_nat)
summary(c1)
#c1_save<-  plot(effect("noncit:race5", c1),   multiline=TRUE, ylab="Avg. Social Reward for Voters", xlab="Proportion Non-U.S. Born Citizens", rug=TRUE, main=NULL)
c1plot <- interact_plot(c1, pred = "noncit", modx = "race5", interval = TRUE, int.width = 0.8, color.class="Set2", x.label = "Proportion Zip Code Foreign-Born", y.label="Social Value") + theme_bw() +
  	 #scale_linetype_manual(values =c(2, 1, 3)) +
	scale_y_continuous(breaks=seq(0.56, 0.84, 0.04)) + expand_limits(y=c(0.56,0.84)) +
	  theme(legend.position = "top",   #Move legend to top
   legend.title=element_blank(), 	# Remove legend name
   legend.key=element_blank(), 	#Remove box around legend points
      legend.text = element_text(size = 28), #increase size of legnd text
   legend.key.width=unit(4,"line"), #increase space between legend symbols
             axis.text=element_text(size=25), #Increase size of axis text
   axis.title.x = element_text(vjust = -0.3, size=28 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5, size=28)	#move x label away from plot
) 
c1plot
 
 pdf('Citizens_Zip_Two_Vote_1plot.pdf', width  = 8, height = 9)
c1plot
dev.off()




#Average Rally Rewards
c2 <- lm(r.mean.stan ~noncit*race5 + inc + edu1 +gender, data=data_LW, weight=wts_nat)
summary(c2)
#c2_save<-   plot(effect("noncit:race5", c2),   multiline=FALSE, ylab="Avg. Social Reward for Ralliers", xlab="Proportion Non-U.S. Born Citizens", rug=TRUE, main=NULL)
c2plot <- interact_plot(c2, pred = "noncit", modx = "race5", interval = TRUE, int.width = 0.8, color.class="Set2", x.label = "Proportion Zip Code Foreign-Born", y.label="Social Value") + theme_bw() +
  	 #scale_linetype_manual(values =c(2, 1, 3)) +
	scale_y_continuous(breaks=seq(0.56, 0.84, 0.04)) + expand_limits(y=c(0.56,0.84)) +
	  theme(legend.position = "top",   #Move legend to top
   legend.title=element_blank(), 	# Remove legend name
   legend.key=element_blank(), 	#Remove box around legend points
      legend.text = element_text(size = 28), #increase size of legnd text
   legend.key.width=unit(4,"line"), #increase space between legend symbols
             axis.text=element_text(size=25), #Increase size of axis text
   axis.title.x = element_text(vjust = -0.3, size=28 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5, size=28)	#move x label away from plot
) 
c2plot

pdf('Citizens_Zip_Two_Rally_1plot.pdf', width  = 8, height = 9)
c2plot
dev.off()
 


#-------- Proportion Latino

#Average Voter Rewards -- Latino vs. White
h1 <- lm(v.mean.stan ~prop.h*race5 + inc + edu1 +gender, data=data_LW, weight=wts_nat)
summary(h1)
#c1_save<-  plot(effect("prop.h:race5", c1),   multiline=TRUE, ylab="Avg. Social Reward for Voters", xlab="Proportion Non-U.S. Born Citizens", rug=TRUE, main=NULL)
h1plot <- interact_plot(h1, pred = "prop.h", modx = "race5", interval = TRUE, int.width = 0.8, color.class="Set2", x.label = "Proportion Zip Code Latino", y.label="Social Value") + theme_bw() +
  	 #scale_linetype_manual(values =c(2, 1, 3)) +
	scale_y_continuous(breaks=seq(0.56, 0.84, 0.04)) + expand_limits(y=c(0.56,0.84)) +
	  theme(legend.position = "top",   #Move legend to top
   legend.title=element_blank(), 	# Remove legend name
   legend.key=element_blank(), 	#Remove box around legend points
      legend.text = element_text(size = 28), #increase size of legnd text
   legend.key.width=unit(4,"line"), #increase space between legend symbols
             axis.text=element_text(size=25), #Increase size of axis text
   axis.title.x = element_text(vjust = -0.3, size=28 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5, size=28)	#move x label away from plot
) 
h1plot
 
 pdf('Latino_Zip_Two_Vote_1plot.pdf', width  = 8, height = 9)
h1plot
dev.off()



#Average Rally Rewards
h2 <- lm(r.mean.stan ~prop.h*race5 + inc + edu1 +gender, data=data_LW, weight=wts_nat)
summary(h2)
#c2_save<-   plot(effect("prop.h:race5", c2),   multiline=FALSE, ylab="Avg. Social Reward for Ralliers", xlab="Proportion Non-U.S. Born Citizens", rug=TRUE, main=NULL)
h2plot <- interact_plot(h2, pred = "prop.h", modx = "race5", interval = TRUE, int.width = 0.8, color.class="Set2", x.label = "Proportion Zip Code Latino", y.label="Social Value") + theme_bw() +
  	 #scale_linetype_manual(values =c(2, 1, 3)) +
	scale_y_continuous(breaks=seq(0.56, 0.84, 0.04)) + expand_limits(y=c(0.56,0.84)) +
	  theme(legend.position = "top",   #Move legend to top
   legend.title=element_blank(), 	# Remove legend name
   legend.key=element_blank(), 	#Remove box around legend points
      legend.text = element_text(size = 28), #increase size of legnd text
   legend.key.width=unit(4,"line"), #increase space between legend symbols
             axis.text=element_text(size=25), #Increase size of axis text
   axis.title.x = element_text(vjust = -0.3, size=28 ), 	#move x label away from plot
     axis.title.y = element_text(vjust = 1.5, size=28)	#move x label away from plot
) 
h2plot

pdf('Latino_Zip_Two_Rally_1plot.pdf', width  = 8, height = 9)
h2plot
dev.off()


 
#Make Table 3 -- for LateX
require(texreg)
texreg(list(b1, b2, h1, h2, c1, c2),
	custom.coef.names = c( 
	"Intercept",
	"Black in Zip code", 
	"White",
	"Income",
	"Education",
	"Male",
	"Zip*White",
	"Citizens in Zip code",
	"Zip*White",
	"Proportion Hispanic",
	"Zip*Hisp"
	),
	custom.model.names = c("Voter", "Rally","Voter", "Rally", "Voter", "Rally" ),
	 dcolumn=TRUE, booktabs = TRUE, use.packages = FALSE, label ="tab:regTable", stars = c(0.05), caption = "Neighborhood Context Effects on Social Value of Partipation", float.pos = "hb", single.row = TRUE) #change significane levels: stars = c(0.01, 0.001)
 

#Add in Means and Standard Deviation for Table Results
require(Hmisc)
#Mean and Variance for Black - Whtie Sample
#Voting
xm <- wtd.mean(data_BW$v.mean.stan, data_BW$wts_nat)
xm
var <- wtd.var(data_BW$v.mean.stan, data_BW$wts_nat)
var
sd <- sqrt(var)
sd
#Rally
xm <- wtd.mean(data_BW$r.mean.stan, data_BW$wts_nat)
xm
var <- wtd.var(data_BW$r.mean.stan, data_BW$wts_nat)
var
sd <- sqrt(var)
sd
#Mean and Variance for Latino - Whtie Sample
#Voting
xm <- wtd.mean(data_LW$v.mean.stan, data_LW$wts_nat)
xm
var <- wtd.var(data_LW$v.mean.stan, data_LW$wts_nat)
var
sd <- sqrt(var)
sd
#Rally
xm <- wtd.mean(data_LW$r.mean.stan, data_LW$wts_nat)
xm
var <- wtd.var(data_LW$r.mean.stan, data_LW$wts_nat)
var
sd <- sqrt(var)
sd





#-- Robustness Check
# Does population density of zip code affect the relationship?
# Does Number of Asian Americans affect outcome?
#Test when addition Asian population control affects the results
summary(c1)
summary(lm(v.mean.stan ~noncit*race5 + inc + edu1 +gender +prop.a, data=data_LW, weight=wts_nat))

summary(c2)
summary(lm(r.mean.stan ~noncit*race5 + inc + edu1 +gender +prop.a, data=data_LW, weight=wts_nat))

summary(b1)
t1<- lm(v.mean.stan ~prop.b*race5 + inc + edu1+gender + UEQE001, data=data_BW, weight=wts_nat)
summary(t1)

summary(b2)
t2<- lm(r.mean.stan ~prop.b*race5 + inc + edu1+gender + UEQE001, data=data_BW, weight=wts_nat)
summary(t2)


summary(c1)
t3<- lm(v.mean.stan ~noncit*race5 + inc + edu1 +gender + UEQE001 + prop.a, data=data_LW, weight=wts_nat)
summary(t3)

summary(c2)
t4<- lm(r.mean.stan ~noncit*race5 + inc + edu1 +gender + UEQE001 + prop.a, data=data_LW, weight=wts_nat)
summary(t4)


summary(h1)
t5<- lm(v.mean.stan ~prop.h*race5 + inc + edu1 +gender+ UEQE001, data=data_LW, weight=wts_nat)
summary(t5)

summary(h2)
t6<- lm(r.mean.stan ~prop.h*race5 + inc + edu1 +gender + UEQE001 ,data=data_LW, weight=wts_nat)
summary(t6)



#Make Appendix Table A3
texreg(list(t2, t1, t6, t5, t4, t3),
	custom.coef.names = c( 
	"Intercept",
	"Black % Zip", 
	"White",
	"Income",
	"Education",
	"Male",
	"Number Residents",
	"Black Zip*White",
	"Latino % Zip",
	"Latino Zip*White",
	"Not US Born % Zip",
	"Not US Zip*White",
	"Prop. Zip Asian"
	),
	custom.model.names = c("Rally", "Voter","Rally", "Voter", "Rally", "Voter" ),
	 dcolumn=TRUE, booktabs = TRUE, use.packages = FALSE, label ="tab:regTable", stars = c(0.05), caption = "Neighborhood Context Results Controlling for Zip Code Population Density", float.pos = "hb", single.row = TRUE) #change significane levels: stars = c(0.01, 0.001)
 





















